THE SELF-GRAVITATING GAS IN THE PRESENCE OF DARK ENERGY: 
MONTE-CARLO SIMULATIONS AND STABILITY ANALYSIS 
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The self-gravitating gas in the presence of a positive cosmological constant A is studied in thermal 
equilibrium by Monte Carlo simulations and by the mean field approach. We find excellent agreement 
between both approaches already for iV = 1000 particles on a volume V [The mean field is exact in 
the infinite N limit]. The domain of stability of the gas is found to increase when the cosmological 
constant increases. The particle density is shown to be an increasing (decreasing) function of the 
distance when the dark energy dominates over self-gravity (and vice- versa). We confirm the validity 
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of the thermodynamic limit: N, V — * oo with N/V^ and A V s fixed. In such dilute limit extensive 
thermodynamic quantities like energy, free energy, entropy turn to be proportional to N. We find 
that the gas is stable till the isothermal compressibility diverges. Beyond this point the gas becomes 
a extremely dense object whose properties are studied by Monte Carlo. 
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I. INTRODUCTION 



The self-gravitating gas in thermal equilibrium has been thoroughly studied since many years0 0, IE S IE S 
OS III EE. As a consequence of the long range attractive Newton force, the selfgravitating gas admits a consistent 
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thermodynamic limit N, V — > oo with fixed. In this limit, extensive thermodynamic quantities like energy, free 

energy, entropy are proportional to iV|^§. 

In ref. Q we investigated how a cosmological constant affects the properties of the non-relativistic self-gravitating 
gas in thermal equilibrium by mean field methods. The mean field approximation becomes exact in the limit when 
the number of particles becomes infinity. 

In the present paper we study the stability properties of the self-gravitating gas in the presence of a cosmological 
constant by mean field and Monte-Carlo methods. The use of Monte-Carlo simulations is particularly useful since 
they are like real experiments and allow to unambiguously determine the stability or instability of the system. 

We compute both by Monte Carlo and mean field methods several physical quantities: the equation of state, the 
particle density and the average particle distance for different values of the dark energy. An excellent agreement 
between both approaches is found except very close to the collapse transition [see figs. 4-6]. The difference between 
both approaches turns to be, as expected, of the order where the number of particles N we choose in the Monte 
Carlo simulations was N = 1000 — 2000. The slightly larger split between Monte Carlo and mean field near the 
collapse is due to the fact that there the corrections to the mean field become singular 0, [j| . 

We find that the onset of instability in the canonical ensemble coincides with the point where the isothermal 
compressibility di verg es. At this point we find that the dimensionless parameter f = G m T i p is maximal. (£ was 
introduced in ref. |l0j for A = 0). We find that the domain of stability of the gas increases for increasing cosmological 
constant. The dark energy has an anti-gravity effect that disfavours the collapse pushing apart the particles. 

In absence of cosmological constant the particle density p(r) is a decreasing function of the distance in the self- 
gravitating gas (see for example @]). In the presence of a positive A we find that p(r) decreases with r when the 
self- gravity dominates over the dark energy. This happens for X = 2A£ < 1- For X > 1 the cosmological constant 
dominates over the self-gravity and the the particle density turns to increase with the distance. This is a consequence 
of the repulsive character of the dark energy. 

The Monte-Carlo study of the self-gravitating gas shows the validity of the dilute thermodynamic limit introduced 
in refs.0,01. Namely, for N, V — > oo with and A V% fixed, a consistent thermodynamic limit is reached where 

extensive thermodynamic quantities like energy, free energy, entropy are shown to be proportional to N. Furthermore, 
the Monte-Carlo simulations allow us to study the condensed phase which turns to be an extremely dense body. The 
phase transition to collapse is found to be of zeroth order. 

The outline of the paper is as follows. Section II presents non-relativistic selfgravitating particles in the presence 
of the cosmological constant, their statistical mechanical treatment and the mean field approach. In section III we 
analyze the stability of the selfgravitating gas for zero and nonzero cosmological constant A while in section IV we 
give the physical picture of the gas and its particle density as a function of A. Section V contains the results of our 
Monte-Carlo simulations on the phase diagram, the collapse transition, the particle density and the properties of the 
collapsed phase. 



II. THE SELF-GRAVITATING GAS IN THE PRESENCE OF THE DARK ENERGY 

A. Non-relativistic selfgravitating particles in the presence of the cosmological constant 

We recall some results about the self-gravitating gas in the presence of the cosmological constant In the 
Newtonian limit the Einstein equations of general relativity become y| 

V 2 V = AirGp-%irGk . (2.1) 

where V stands for the gravitational potential, p for the density of massive particles and A for the cosmological 
constant. Eci. (|2.1|) determines the weak and static gravitational field produced by non relativistic matter in the 
presence of the cosmological constant. For zero cosmological constant we recover the usual Poisson equation, as it 
must be. 

The Hamiltonian for such set of self-gravitating particles in the presence of the cosmological constant can be written 
in the center of mass frame as£| 

where rrn stands for the mass of the particle at the point qi with momentum pi and 

2J rrij ql = , 
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since we choose center of mass coordinates. 

The cosmological constant contribution to the potential energy grows negative for increasing values of the particle 
distances qi to the center of mass. Therefore, the gravitational effect of the cosmological constant on particles amount 
to push them outwards. It can be then said that the cosmological constant has an anti-gravitational effect. 



B. Statistical mechanics of the self-gravitating gas with A ^ 

We present here the statistical mechanics of the self-gravitating gas in the presence of the cosmological constant in 
the canonical ensemble. For simplicity we shall consider TV particles with identical mass m. 

The Hamiltonian is given by ea. I|2.2|l and therefore the classical partition function of the gas is then, 

It is convenient to introduce the dimensionless coordinates ^ = 77 ■ The momenta integrals are computed straightfor- 
wardly. Hence, the partition function factorizes as the partition function of a perfect gas times a coordinate integral. 



z =^(£) l-JUfn^^, (2.3) 



where 

N 



l<i<j<N 1 1 31 i=l 

The dimensionless parameters which characterize the strength of self-gravity and dark energy are respectively, 

Gm 2 N L 2 , . 

r ] =^ F1 - and £ = 2AGm — . (2.5) 

77 can be interpreted as the ratio of the self-gravitating energy per particle to the kinetic energy per particle. 
We define the ratio of dark energy and self-gravity 

77 ml\ 

It is the ratio between the dark energy contained inside the volume V and the mass of ordinary matter. 

In the large N limit the 3iV-uple integral in eq. (|2.3|l can be approximated by a functional integral over the density. 
In the saddle point approximation we find that the density 



obeys to the differential equation £| 



The density has to be normalized as 



P (x) = e*W . (2.7) 
V 2 $(x) + 4tt ( 77 e*« - t) = . (2.8) 



J d 3 xp(x) = 1 . (2.9) 



In hydrostatic equilibrium we obtain the same equation 0. Therefore, hydrostatics and mean field are equivalent in 
the N — > 00 limit. 

In the limiting case X = ea . (12.81) b ecomes the well known Lane-Emden equation in the absence of cosmological 



in tne limiting case A — u ca. iiz.»|) becomes tl 
constant (see refs. H H H S H EMI 13 ) 



4 



1 r 

f vs. eta for X = 

f vs. eta for X = 0.3 - 

f vs. eta for X = 1 
f vs. eta for X = 1 .5 



1.4 



1.2 - 



0.2 




0.5 



1.5 



2.5 



eta 



FIG. 1: The density at the boundary f(X,rj) versus r\ for X = 0,0.3,1,1.5 by the mean field approach. The instability 
points ri (X) where the thermal compressibility diverges are pinpointed on the plot by a +. They are r\ (X = 0) = 1.510 . . ., 
r) (X = 0.3) = 1.63 . . ., r/ (X = 1) = 2.04 ... and r] (X = 1.5) = 2.55 . . .. The gas is stable from t] = till r] = r) (X) in the 
canonical ensemble. 



C. The isothermal sphere with A / 

For spherically symmetric configurations the mean field equation l|2.8[) becomes an ordinary non-linear differential 
equation. 



d 2 $ 2 d$ 



dR 2 



(2.10) 



r-) 3 . The density of particles p(R) has to be 



The various thermodynamic quantities are expressed in terms of their solutions. We work in a unit sphere, therefore 
the radial variable runs in the interval < R < R m ax, Rmax 

normalized according to eq. (|2.9() . Integrating eq. I|2.10|l from R = to R — R max yields, 

ri-Z = -R 2 max $'(JU) 

Setting, 



(2.11) 



$(i?) = u(x) + In ^-s 



'3? 



R 



Rr, 



(2.12) 
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£ R = jj£ — and T] R — — , the saddle-point equation (|2.10|l simplifies as, 
In order to have a regular solution at origin we have to impose 



d 2 u 2 du 
dx 2 x dx 



u'(0) = . (2.14) 
We find from eas. (|2.11|) - (|2.12|) for fixed values of r\ and £, 

v R - t R 



'(V^)=-^- (2.15) 



/3£ R 

Eas. H2.14|) and l|2.15|l provide the boundary conditions for the nonlinear ordinary differential equation l|2.13|l . In 

particular, they impose the dependence of uq = u(0) on rj R and £ R . 

We find from ea. l|2.13|l for small x the expansion of u(x) in powers of x with the result, 

u{x) = u + (1 - e u °) y + e"°(e"° - 1) ^ + 0(a; 6 ) . (2.16) 
where we imposed ea. (|2.14() . The value of uo follows by imposing eq. l|2.15[) . 

Using eas. (!2.6|l . I|2.7|l and [|2.12JI we can express the density of particles in terms of the solution of ea. l|2.13|l 



p(R) — X exp 



R 



(2.17) 



We choose to express the physical quantities in terms of rj [ea. (|2.5|l ] and the ratio X of dark energy to self-gravity 
[ea. p.Gfl ]. The local pressure obeys the ideal gas equation of state in local form P(R) = p(R) . Using eas. (|2.17() 
the pressure at the boundary is given by 



NT 



V 



N T 



P{R max ) = -ft- X exp u{ ^/3Xr, R ) = — f(X, V ). (2.18) 



V 



The contrast C(X, rf) between the pressure at the center and at the boundary can be written as, 

i \-tl-max ) 

We plot in fig. ^the density at the boundary versus rj for different values of X. 

III. STABILITY OF THE SELF-GRAVITATING GAS 
A. Stability for A = 0. 

The self-gravitating system can be in one of two phases: gaseous or highly condensed. Mean field theory describes the 
gas phase and gives exactly the physical quantities in the thermodynamic limit. The equilibrium density p{x) solution 
of the mean field equation l|2.8|l minimizes the Helmholtz free energy in the canonical ensemble. This density is the 
most probable distribution which becomes absolutely certain in the infinite N limit. All thermodynamics quantities 
follows from this density p(x). This continuous density p(x) solution of the mean field equation l|2.8|l fails to describe 
the condensed phase. 

The condensed phase have been found by Monte-Carlo simulations for A = [2j. In the microcanonical ensemble 
the center collapses in a very dense core surrounded by a halo of particles. In the canonical ensemble all the particles 
collapse in a very dense body . 

The stability of the gaseous phase in the mean field approach can be analyzed looking for the extrema of the 
dimensionless parameter 

G 3 rn 6 N 2 P 

C= y f ■ (3-1) 



6 




FIG. 2: The value of the parameter rj at the Jeans instability: r\ — rj (X) versus X. 



when the temperature T, the pressure P and the number of particles N are kept fixed in the canonical ensemble. The 
gas phase becomes instable (for A = 0) when £ is maximal jlO|. 



Using eqs.JUHJ, (|^TH|> and (|3~T|) we find that 



where fo(r]) = f(X — 0, 77) is the external density in absence of dark energy. 

A better physical insight is obtained by looking to the the isothermal compressibility 



K T = 



1 fdV 



V \dP 



(3.2) 



(3.3) 



Using eqs. H2.5JI and H2.18|l we find that the dimensionless isothermal compressibility kt = ^r- Kt takes the form 

1 



We see from eqs.J221> an d (|3~I|) that 



(3.4) 



dr/ 



7 



Hence, kt diverges at the extrema of £. 

The isothermal compressibility is positive from 77 = till 77 =j]o = 1-510. . .. At this point kj- as well as the 
specific heat at constant pressure diverge and change their signs0,0|- Moreover, at this point the speed of sound 
at the center of the sphere becomes imaginary 3 . Therefore, small density fluctuations will grow exponentially in 
time instead of exhibiting oscillatory propagation. Such a behaviour leads to the Jeans instability and collapse 14] . 
Monte-Carlo simulations confirmed the presence of this instability at r\ — r\ = 1.510 ... in the canonical ensemble for 
A = 2]. 



B. Jeans instability for A 7^ 



We now compute the dimensionless isothermal compressibility for the self-gravitating gas with A ^ 0. Introducing 
the dimensionless parameter 

2 G 3 m 5 A 2 A 3 

a = — = rr X, (3.5) 

the dimensionless isothermal compressibility kt takes the form: 

k t = 7 — r • (3.6) 

f( V ,a) + § (fj) (*?,«) 

The Jeans instability happens when the isothermal compressibility diverges. This happens for 77 = rj {X) verifying: 

f(Vo, a) + ^ (^-] (Vo, a) = . 



3 \dr) 

Thus, for a fixed a, r\ is the value of 77 which maximizes the quantity 

C = 77 3 /(77,a) = — 4 . (3.7) 

As for the self-gravitating gas with A = 0, C is maximal when the specific heat at constant pressure diverges and 
changes its sign. We plot ?7 versus X in fig. fusing ea. (|3.5jl . Monte Carlo simulations confirm that the onset of the 
Jeans collapse happens at rj = r/ . 
We can write 77 [ea. l|2.5|) ] as 

_ GmM 
V ~ TL 

where M is the total mass of the gas. Therefore, the collapse condition 77 > 770(A) can be written as 

LT 



M > M rjo(X) where M 



G m 



The gas collapses when its mass takes a larger value than the critical one Mj(X) = M rj (X). We can consider 
Mj(X) as the generalization of the Jeans mass in the presence of the cosmological constant. Notice that the presence 
of A increases the Jeans mass with respect the A = case [see fig. 



IV. PHYSICAL PICTURE 



A. Behaviour of the gas phase with X 



The effects of self-gravitation and dark energy go in opposite directions. Self-gravitating forces are attractive, 
while dark energy produces repulsion. If the ratio of dark energy and self-gravity is X < 1 [defined by ea. (|2.6l) ]. 
self-gravitation dominates over dark energy. If X = 1 the effect of self-gravitation is exactly compensated by the dark 
energy. If X > 1 dark energy dominates over self-gravitation. We illustrate these three cases plotting /(A, 77) versus 
t] for X = 0.3, 1, 1.5 in fig. □ 
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• first case: X < 1: self-gravity dominates 

We illustrate this case choosing X = 0.3. We see in fig. ^that /(X, 77) versus 77 exhibits a form analogous to 
f(X — 0,77) versus 77 in absence of dark energy. f(X,rj) has two Riemann sheets as a function of 77. In the first 
sheet f(X,i]) monotonically decreases with 77 for fixed X < 1. 

• second case: X = 1: exact compensation 

The effect of self-gravitation is compensated by the effect of dark energy. We see in fig. Q]that f(X = 1, 77) = 1. 
An exactly homogeneous sphere p(R) = 1 is a solution (see below). In such special case the self-gravitating gas 
behaves as a perfect gas (with PV = NT everywhere). 

• third case: X > 1: dark energy dominates 

We illustrate this case choosing X = 1.5. We see in fig. [flthat f(X = 1.5, 17) is an increasing function of r\ for 
fixed X > 1. 

The Jeans instability happens in the three cases. The value of 77 = r/ (X) where the gas collapses increases with X 

(fig. The gas collapses for X > at a higher density for a given temperature or for a lower temperature for 

v 3 

a given density than for X = 0. That is, the domain of stability of the gas increases for increasing X . The dark 

v 3 

energy has an anti-gravity effect that disfavours the collapse pushing the particles toward the boundary of the sphere. 
Notice that the point of Jeans instability is in the first Riemann sheet of the 77-plane for < X < 1. 



B. The particle density p(R) 

We study the density p(R) as a function of the radial coordinate R. The density p(R) of the self-gravitating gas in 
absence of dark energy (X = 0) is always a decreasing function of R. This follows from the attractive character of 
gravitation (see for example ref. 0,0). 

The i?-dependence of the density p{R) of a self-gravitating gas in the presence of dark energy is more involved 
because the dark energy opposes to the attraction by gravity. The functions p{R) and u(x) have the same qualitative 
dependence in R since they are related by exponentiation [see eq. I|2.17|l ] . 

The behaviour at the center of the sphere (R = 0) is governed by the sign of tj since u'(0) = [eq. (|2.14[l ] and we 
find from eq.lEHJ that u"(0) = §(1 - e u °). Hence, 

sign[u"(0)] = -sign[u ] . 

Therefore, p(R) decreases (increases) at R — if u > (u < 0). 

The behaviour at the boundary of the sphere (R = R m ax) is governed by the sign of 7/ — £ since according to 
ea. (l2~T5l) 

sign[u'(R max )} = -siga[r) R - £ R ] = sign(X - 1) . 

Therefore, p(R) decreases (increases) at R ma x if X < 1 (X > 1). 

The particle density exhibits for the stable solutions (namely, for kt > 0) one of the following behaviours illustrated 
in fig. 

• decreasing: X < 1 and uq > 0. The density p(R) decreases from the center of the sphere till the boundary. 
We plot in fig. [21 the density p(R) versus R for X — 0.3 and rj = 1. The self-gravitation dominates over dark 
energy. 

• increasing: X > 1 and uq < 0. The density p{R) increases from the center of the sphere till the boundary. We 
plot in fig. [3]the density p(R) versus R for X = 1.5 and 77 = 1. The dark energy dominates over self-gravitation. 

• homogeneous: X = 1 and uq — 0. The density is homogeneous: p(R) = 1. The effect of self-gravity of 
particles is exactly compensated by the effect of dark energy. We plot in fig. |3]the density p(R) versus R for 
X = 1 and 77 = 1. 

Eq. 12.8fl for X = 1 is similar to the equation that describes in two space dimensions (multi)-vortices in the Ginsburg- 
Landau or Higgs model in the limit between superconductivity of type I and II |17|. However, for the vortex case one 
has £ = 77 < since like charges repel each other in electrodynamics while masses attract each other in gravity. This 
opposite nature of the forces is responsible of the fact that eq. I|2.8|l for X = 1 only has trivial stable solutions while a 
host of non- trivial solutions appear in the vortex case|l7j. 
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rho(R) vs. R for X=0.3 and eta=1 
rho(R) vs. R for X=1 and eta=1 
rho{R) vs. R for X=1 .5 and eta=1 
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FIG. 3: The particle densities p(R) vs. the radial coordinate R for X = 0.3 and r\ = 1 (decreasing), X — 1 and 77 = 1 
(homogeneous) and finally X — 1.5 and rj = 1 (increasing) computed by the mean field approach 

V. MONTE CARLO CALCULATIONS 
A. The Metropolis algorithm 

The standard Metropolis algorithm 0] was first applied to self-gravitating gas in absence of dark energy in ref.0. 
We apply this standard Metropolis to the self-gravitating gas in presence of dark energy. We perform it in a volume 
V in the canonical ensemble at temperature T. We compute in this way the external pressure, the energy, the average 
density, the average particle distance and the average squared particle distance as function of r] for a given value of 
the ratio of dark energy and self-gravity X = | = • 

We implement the Metropolis algorithm in the following way. We start from a random distribution of N particles in 
the chosen volume. We update such configuration choosing a particle at random and changing at random its position. 
We then compare the energies of the former and the new configurations. We use the standard Metropolis test to 
choose between the new and the former configurations. The energy of the configurations are calculated performing 
the exact sums as in ea. H2.4L We use as statistical weight for the Metropolis algorithm in the canonical ensemble, 



[see eas. (|2.3|) and 1)2.4)) ] . The number of particles go up to 2000. 
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We introduce a short distance cutoff in the Newtons' potential [see ea. (|2.4|l ] 

. 1 1 



l<i<]<N 11 



with 



\n-fj\A = \n-fj\ , for \j f i -fj\>A 
= A , for \f t - fj\ < A , 

where A <^i_V^ is the short distance cut-off. The presence of the short distance cut-off prevents the collapse (here 
unphysical) of the self-gravitating gas (for more details see 0). 

B. Influence of the geometry 

We apply the Metropolis algorithm to the self-gravitating gas in presence of dark energy in a sphere and in a cube. 
For X = (absence of dark energy) the Monte Carlo results in a sphere and in a cube with identical volume give 
practically the same results |2| ■ The geometry does not influence the physics of the system. 

However, for X > (presence of dark energy) the Monte Carlo results in a sphere and in a cube turn to be different. 
This effect can be traced to the repulsive character of the dark energy that pushes the particles towards the boundary. 
As a result, the physical quantities becomes sensible to the geometry. We thus compared the mean field results with 
spherical symmetry from the previous sections and ref.[4| with Monte Carlo calculations performed on a sphere. 



C. Phase diagrams 

We recall that our self-gravitating system exhibits three different behaviours according the value of X, the relative 
rate of dark energy and self-gravity. In all three cases two different phases show up: for r\ < t]t(X) we have a non 
perfect gas and for 7/ > t]t(X) it is a very condensed system with negative pressure. The transition between the two 
phases is very sharp. 

• For X < 1 the self-gravity dominates over dark energy. We plot the external density f(X, rf) = versus r\ 
for X = 0.3 in fig. In the gas phase monotonically decreases with rj from rj = till the collapse point 
r] = T] T (X). 

The average distance between particles < \fi — Tj\ > and the average squared distance between particles 
i/< \ fi — rj\ 2 > monotonically decrease with i] (fig. When the gas collapses at jj^, < \fi — fj\ > and 

< \ fi — fj\ 2 > exhibit a sharp drop. 

• For X = 1 the effect of dark energy exactly compensates the effect of self-gravity. We have — 1 in the gas 
phase. The system behaves like a perfect gas from rj = till the collapse point rj — i]t(X). 

• For X > 1 the dark energy dominates over self-gravity. We plot in fig. |S]the external density f{X,rj) = 
versus i] for X — 1.5. In the gas phase monotonically increases with rj from r\ = till the collapse point 

T] = T] T (X). 

We find that the Monte Carlo calculations in the sphere accurately reproduce the mean field results for in the 
gas phase with spherical symmetry. We then consider a cube and a sphere of identical volume. The Monte Carlo 
calculations turn to give different results in these two cases: a) The external pressure is larger in the cube than in the 
sphere, b) The gas phase is more stable in the cube, c) The collapse value t]t{X) in the cube is larger than its value 
for the same X in the sphere. 



D. The collapse point 

The phase transition happens in the Monte Carlo simulations at r/ = i]t(X). It must be compared with the point 
of Jeans instability rj = r) (X) according to mean field [see III.C] where the isothermal compressibility diverges and 
changes its sign. We see in that t]t(X) is very close of r/ (X). They are probably the same point. For X = 0.3 
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FIG. 4: The density at the boundary fx (if) = versus 77 for X = 0.3 by Monte Carlo methods and mean field approach. 
The collapse value of the gas phase is rfr(X = 0.3) = 1.63 . . . 

we have Tj (X) = 1.63 ... and T] T (X) = 1.63 .... For X = 1.5 we have rj a (X) = 2.55 ... and T) T (X) = 2.6 . . .. We 
conclude that the collapse of the system observed in the Monte Carlo simulations is due to the Jeans instability when 
the isothermal compressibility diverges and is well described by the mean field approach. 

For 77 > t)t(X) we see that the gaseous phase is metastable. Before collapsing the system stay for a long Monte 
Carlo time in the metastable gaseous phase in the Monte Carlo computations. The Monte Carlo time for collapse 
increases with increasing value of X. This is due to the repulsive effect of dark energy. The gas phase is less unstable 
in the presence of dark energy. 

E. Average distribution of particles 

We illustrate the two behaviours of the self-gravitating gas in the presence of the cosmological constant plotting 
the average density as a function of two Cartesian coordinates, the third coordinate being fixed. For simplicity we 
depict the densities in the cube. 

For X > 1 we see that the density is larger on the boundary that at the center of the cube (fig. 0) . The dark energy 
dominates over the self-gravity and pushes the particles toward the boundary of the cube. 

For X < 1 we see that the density is larger at the center than on the boundary of the cube (fig. [SJ. The self-gravity 
dominates over dark energy and attracts the particles to the center of the cube. 

We illustrate now the two phases (gaseous and condensed) plotting the average particle distribution in a cube. Fig. 
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FIG. 5: The average particle distance and the average squared particle distance versus r\ for X = 0.3 by Monte Carlo methods. 
They exhibit a sharp decrease at the collapse value r\ = t]t{X — 0.3) = 1.63 . . .. 

Eland fig. EH depict the average particle distribution from Monte Carlo calculations with 1000 particles for X = 0.3 
at both sides of the collapse point. Fig. corresponds to the gaseous phase and fig. ^|to the collapsed phase. 

We find that the Monte Carlo simulations (describing thermal equilibrium) are much more efficient than the TV-body 
simulations integrating Newton's equations of motion. [Indeed, the integration of Newton's equations provides much 
more information than thermal equilibrium investigations] . Actually, a few hundreds of particles are enough to get 
quite accurate results in the Monte Carlo simulations (except near the collapse points). Moreover, the Monte Carlo 
results turns to be in excellent agreement with the mean field calculations up to corrections of the order 



F. The condensed phase 

In the condensed phase all the particles collapse in a very dense body (see fig. I10[l . The self-gravity contribution 
of the potential energy dominates overwhelmingly the dark energy contribution in the condensed phase. Using eqs. 
(|2.2|) . H2.5|) and l|2.6|) and the virial theorem Q we obtain that the external pressure is expressed as 



/x(r?) = 



PV 
NT 
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PV/NT vs. eta for X=1 .5 Monte Carlo method 
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FIG. 6: The density at the boundary fx(v) = versus tj for X = 1.5 by Monte Carlo methods and mean field approach. 
The collapse value of the gas phase is rfr(X = 1.5; = 2.6 . . . 



The Monte Carlo results indicate that 
approximated by 



50. Thus in this condensed phase the external pressure can be 



PV 



where K ~ 16. 

Since f(rf) has a jump at the transition, the Gibbs free energy is discontinuous and we have a phase transition of 
the zeroth order as in the absence of cosmological constant Q ■ 



VI. DISCUSSION AND CONCLUSIONS 



The behaviour of the self-gravitating gas is significantly influenced by the cosmological constant (or not) depending 
on the value of the ratio X defined by ea. l|2.6|l . 

2A7 / cosmological constant \ A 
A = = 2 I I = 2 — 

m N \ mass / v py 

where py = N m/V is the average mass density in the volume V and where the subscript () v indicates mass and 
cosmological constant inside the volume V. 
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X takes small values for astrophysical objects except for clusters of galaxies where X galaxy cluster s ~ 0.01 and 
especially for the universe as a whole where Xuniverse ~ 4. The smallness of X stems from the tiny value of the 
present cosmological constant A = 0.663 10 _29 g/cm 3 18]. The particle mass density is much larger than A in all 
situations except for the universe as a whole. 

However, the non-relativistic and equilibrium treatment does not apply for the universe as a whole. Certainly, 
the galaxy distribution can be considered as a non-relativistic self-gravitating gas in the presence of the cosmological 
constant. But the gas of galaxies is certainly not yet at thermal equilibrium today. Perhaps, thermal equilibrium 
may be reached in some billions of years. Incidentally, at thermal equilibrium the particle density is regular at the 
center of the distribution in the whole range of X. Namely, there is no cusp at the center in a situation of thermal 
equilibrium. 




FIG. 7: The local density of the gas p(x, y, z) versus the Cartesian coordinates x and y at z = 5 in a cube of size 9 for X — 3 
and rj — 1 from Monte Carlo simulations. The density is larger on the boundary of the cube than at the center. 
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rho(x,y,z) vs. x and y at z=5 for X=0.3 and eta=1 




FIG. 8: The local density of the gas p(x,y,z) versus the Cartesian coordinates x and y at z — 5 in a cube of size 9 for X = 0.3 
and ij = 1 from Monte Carlo simulations. The density is larger at the center of the cube than on the boundary. 
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FIG. 9: Average particle distribution in the gaseous phase from Monte Carlo simulations in a cubic volume for X = 0.3, r/ = 1.6 
and N = 2000. A + denotes one particle. 
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FIG. 10: Average particle distribution in the collapsed phase from Monte Carlo simulations in a cubic volume for X = 0.3, 
r\ = 1.8 and N = 2000. A + denotes one particle. 



